% 原始代码

clear
clc
N = 4; M = 4;
H= (randn(M,N) +j*randn(M,N))/sqrt(2);
Pathloss = 1; 
H = Pathloss*H;
SNR_dB = -10:5:30;

for I_SNR = 1 :length(SNR_dB)
    rho = 10^(SNR_dB(I_SNR)/10);
    F = rho*H'*H;
    [U_F,L_F] = eig(F); % 用于计算矩阵的特征值和特征向量,U_F是特征向量矩阵，L_F是特征值对角矩阵。
    
    up_u = 10^10;  lo_u = 0;   err = 10;
    Lambda_F = max(real(diag(L_F)),0);
    while ( abs(err) > 10^(-5))
        u = (up_u + lo_u)/2;
        p_p = max(u - 1./Lambda_F,0);
        err = abs(trace(U_F*diag(p_p)*U_F')) - 1;
        if (err > 0)
            up_u = u;
        else
            lo_u = u;
        end
    end
    optQ = U_F*diag(p_p)*U_F';
    
%     for循环为计算均值，循环1000次，然后求平均   
%     for t = 1:1000
%     Capacity_temp(t) = log2(abs(1+rho*h(t)*h(t)') );
%     end
%     Capacity(I_SNR)= sum(Capacity_temp)/NUM_timesamples;

    Capacity(I_SNR) = log2(abs(det(eye(M) +rho*H*optQ*H') ));
    Rate_EP(I_SNR) = log2(abs(det(eye(M)+rho*H*eye(N)/N*H')) ) ;
    
    [value_LF,Id_LF] = max(diag(L_F)); 
    p_bf = zeros(N,1);
    p_bf(Id_LF) = 1;   BF_Q = U_F*diag(p_bf)*U_F ;
    Rate_BF(I_SNR) = log2(abs(det(eye(M)+rho*H*BF_Q*H')) );
end

figure(1)
hold on
plot(SNR_dB, Capacity, 'ro-' , SNR_dB, Rate_EP, ' b^-',SNR_dB,Rate_BF,'k*-')
legend('Optimal scheme' , 'Equal power allocation' ,'Beamforming' ,'Location' , ' NorthWest'),
xlabel(' rho (in dB)')
ylabel('Capacity (bps/Hz)')
grid

